Transcriptomes reveal the involved genes in the sea urchin Mesocentrotus nudus exposed to high flow velocities

Despite the importance of flow velocity in marine ecosystems, molecular mechanisms of the water flow induced behavioral and growth changes remain largely unknown in sea urchins. The present study compared the gene expressions of the sea urchin Mesocentrotus nudus at high flow velocities (10 cm/s and 20 cm/s) and low flow velocity (2 cm/s) using transcriptomes. A total of 490 and 470 differentially expressed genes (DEGs) were discovered at 10 cm/s and 20 cm/s, respectively. There were 235 up-regulated and 255 down-regulated genes at 10 cm/s, 213 up-regulated and 257 down-regulated genes at 20 cm/s, compared with sea urchins at 2 cm/s. Further, there were 72 overlapped DEGs involved in regulation at both 10 cm/s and 20 cm/s. Gene Ontology (GO) functional annotation showed that DEGs were mainly enriched to cellular process, cell part, binding, and metabolism process. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis found that DEGs were enriched in three pathways related to amino acid metabolism and lipid metabolism. A number of genes related to growth and metabolism of sea urchins were mobilized in high flow velocity environment. We further highlighted a muscle-associated gene ankyrin-1, which is correlated with the movement of tube feet at different flow velocities. The present study provides valuable information on the molecular mechanisms of changed behaviors and growth when sea urchins are exposed to high flow velocity.


Results
Transcriptome sequencing and assembly. A total of 432,777,420 raw reads were sequenced from nine samples. 403,298,480 (93.17%) clean reads were assembled into 156,177 unigenes, in which the average length was 796 bp, mainly distributed in 200-1200 bp (Fig. 1A) and GC-content was 38.39%. There were 82.39% of the reads (271,194,264) compared to the transcript sequences. Benchmarking Universal Single-Copy Ortholog (BUSCO) (v3.0.1) evaluation was performed on the assembled unigenes and 94.72% of transcripts were compared with BUSCO groups.

Gene function annotation and classification. The comparison results of unigenes in various databases
are shown in Table 1. There were 37,250, 18,934, and 24,546 unigenes annotated Nr (NCBI non-redundant protein sequences), GO (Gene Ontology), and Swiss-Prot (a manually annotated and reviewed protein sequence database), respectively (Fig. 1B).
The transcriptome assembly was annotated and categorized based on GO (Fig. 1C). In the biological process category, unigenes were mainly annotated in cellular process, metabolic process, and biological regulation. In cell composition, they were mainly concentrated in cell part and organelle. The molecular functions were mainly enriched in binding and catalytic activity.
KOG annotation results are showed in Fig. 2A. Most genes were annotated in general function prediction only, signal transduction mechanisms, and posttranslational modification, protein turnover, associated protein. Real-time quantitative PCR validation of RNA-Seq data. Real-time quantitative PCR and transcriptome sequencing of differential expression multiple results are shown in Fig. 5. In the validation of 12 DEGs, qRT-PCR results showed the same change trend as the transcriptome results, indicating that the transcriptome results were reliable.

Discussion
Mesocentrotus nudus mainly inhabit in shallow seas and intertidal zones with large changes in water flow 27 . Growth and movement commonly show great differences in M. nudus at different flow velocities 28 . However, it remains largely unclear on the molecular mechanism of the effects of flow velocities on those performances of M. nudus. The present study observed 403,298,480 clean reads and Q30 over 93.64% by transcriptome sequencing. A total of 156,177 unigenes were obtained from tube feet, which were less than that in tube feet of the sea urchin Strongylocentrotus intermedius 29 , but more than that from gonads of M. nudus 2 . This suggests that the present sequencing results are qualified and provide valuable resources for the research on the functional genes of M. nudus. The gene expression analysis found that sea urchins exposed to high flow velocity (20 cm/s) showed 470 DEGs, which are obviously less than in other environmental stressors. For example, 2125 DEGs were identified in gonads of S. intermedius at high water temperature 30 and 29,107 DEGs were found in coelomocytes of S. intermedius in a hypoxia environment 31 . Different environmental stressors induce constant expression regulation of functional genes and consequently contribute to the cellular homeostasis 32 . The number of DEGs was less when sea urchins were exposed to a high flow velocity, indicating that high flow velocity probably causes a lower stress to sea urchins than previously being thought. The regulatory mechanisms of M. nudus may be relatively simple in response to high flow velocity. Notably, the gonad is a tissue that possess more complex metabolism than tube feet. More uniquely expressed genes were found in gonads than in tube feet of the sea urchin Tripneustes gratilla 33 . Thus, the complexity of the regulatory mechanisms is not exclusive to explain the large difference in the number of DEGs and the main reason needs to be further studied. Further, there are only 72 overlapped DEGs involved in regulation at both 10 cm/s and 20 cm/s, compared with the total number of DEGs between the groups of 10 cm/s (490 DEGs) and 20 cm/s (470 DEGs). They are important resources for further investigation.  www.nature.com/scientificreports/ There were 235 up-regulated and 255 down-regulated genes at 10 cm/s, and 213 up-regulated and 257 downregulated genes at 20 cm/s, compared to at 2 cm/s. GO annotation showed that the DEGs of sea urchins were enriched in metabolic process, biological regulation and development process. Organisms make metabolic adjustments in response to different environmental stressors 30 . High flow velocity significantly affects growth and metabolism of fish. For example, it greatly increases the energy consumption of the fish Rhynchocypris lagowskii 34 and reduces the fat content of the flounder Paralichthys olivaceus 35 . Consistently, our previous studies found that high flow velocity displays negative effects on the growth and development of sea urchins 18,28 . Here, we further found three metabolism pathways through KEGG analysis, including arachidonic acid metabolism, linoleic acid metabolism, and arginine and proline metabolism. Arachidonic acid metabolism and linoleic acid metabolism are important pathways of lipid metabolism 34 . The arginine synthesis pathway regulates the lipid metabolism through promoting oxidative decomposition of glucose and fatty acids and lipolysis of fat cells 36 . These results www.nature.com/scientificreports/ suggest that pathways related to lipid anabolism probably play an important role in regulating growth and development of sea urchins at high flow velocity. Adhesion and swing of tube feet are the basis to keep balance in the foraging of sea urchins. Time for adhesion significantly extends in sea urchins at high flow velocity, because water flow reduces the swing of tube feet 28 . In the present study, 420 DEGs were found by comparing the genes expression between 10 cm/s and 20 cm/s. Transcriptome data further showed that the expression of ankyrin-1 was significantly higher at 20 cm/s than that at 10 cm/s. The expression of nectin-1 related to adhesion showed significant difference in sea urchins at different flow velocities, which is considered to be the key gene reducing the movement capability of sea urchins at high flow velocity 23 . Ankyrin-1 is annotated to the sarcoplasmic reticulum, which is a master regulator of muscle contraction 37 . Ankyrin-1 was firstly found in the preparation of erythrocyte membranes, which is an erythrocyte membrane protein that connects the underlying cytoskeleton to the plasma membrane 38 . It is consistently expressed in muscle and is annotated in other muscle activity-related GO classifications, including the Z-line 20,39 . Therefore, it is reasonable to speculate that ankyrin-1 probably plays an important role in the movement capacity of sea urchins at different flow velocities. Higher expression of ankyrin-1 is probably involved in mobilizing more energy of muscle to resist the stress of high flow velocity in collaboration with other genes, for example nectin-1. This novel finding increases our understanding on how sea urchins adapt to high flow velocity.

Conclusion
The present study investigated the gene expressions of sea urchins at different flow velocities using transcriptomes. There were 490 and 470 DEGs at flow velocity of 10 cm/s and 20 cm/s compared to 2 cm/s, respectively. GO annotation found that a number of DEGs enriched in growth and metabolism. KEGG pathway analysis discovered three pathways related to lipid anabolism and amino acid metabolism. Thus, tube feet show different gene expression levels when sea urchins live in the habitats at high flow velocity. In addition, we found that 72 overlapped DEGs involved in regulation at both 10 cm/s and 20 cm/s. Functions of these 72 DEGs are worth further investigating to find out the molecular mechanisms of the adaptation of sea urchins exposed to high flow velocity. We highlighted a muscle-associated gene ankyrin-1, which is correlated with the movements of tube feet at different flow velocities. The present study provides valuable information into the molecular mechanisms of changed behaviors and growth of sea urchins exposed to high flow velocities.

Materials and methods
Sea urchins. Mesocentrotus nudus of 8-month age (test diameter of 21.57 ± 0.14 mm, test height of 10.12 ± 0.10 mm and body weight of 4.19 ± 0.08 g) were transported from a local aqua-farm to the Key Laboratory of Mariculture & Stock Enhancement in North China′s Sea, Ministry of Agriculture and Rural Affairs, Dalian Ocean University (121° 56′ E, 38° 87′ N). All sea urchins were then maintained in fiberglass tanks (~ 1000 m 3 ) with the natural photoperiod with aeration and fed fresh kelp Saccharina japonica ad libitum. One third of seawater was renewed and feces were removed every 2 days. The water quality was measured regularly using a water quality monitoring meter (YSI, Incorporated, OH, USA). The salinity was 31.96 ± 0.05‰ and the water temperature was 19.15 ± 0.85 °C.

Experimental design.
A circular runway (Fig. 6) with five rooms as experimental areas (length × width × height = 20 × 15 × 10 cm) was used in the present study 28 . The ball valve and flow meter (JDC Electronic SA Co., Switzerland) were used to control and measure flow velocity. Water circulation was supported by a water pump (200 W, 20,000 L h −1 , Jebao, China). Sea urchins (test diameter of 21.57 ± 0.14 mm) were divided equally into three groups and cultured at three different flow velocities (2, 10, and 20 cm/s) from 7 May 2019 to 17 July 2019. There were 25 sea urchins at each flow velocity, with 5 individuals in each of the 5 rooms.
Six sea urchins were randomly selected from all sea urchins in the device. The tube feet from two sea urchins were collected in a 1.5 mL centrifuge tube and mixed as one sample to meet the test standard (12 µg) at the end of the experiment. There were three biological replicates for each experimental group (N = 3). The nine samples were immediately frozen in liquid nitrogen and stored at − 80 °C. RNA extraction, library construction, and RNA sequencing. Total RNAs were extracted using the conventional phenol/chloroform extraction method 40 . RNAs were diluted in a certain proportion. A microspectrophotometer (K5500, Beijing) and the Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA) were subsequently used to measure RNA concentration, purity, and integrity, respectively. The sample standard was OD260/280 ≥ 1.8 and OD260/230 ≥ 1.8.
The qualified samples were used for RNA library construction 41 . The mRNA was enriched by Oligo (dT) beads. The fragmented mRNA was used as the template to produce the first-strand cDNA. The first-strand cDNA was added in fragmentation buffer, dNTPs, RNaseH, and DNA Polymerase I to produce the second-strand cDNA. The products were used for the size selection by gel electrophoresis after being purified using QiaQuick PCR Kit (Qiagen, Germany). Finally, the products were amplified by PCR. The library was sequenced using Illumina Novaseq 6000 with the sequencing strategy of PE150.
The raw reads were processed through Perl scripts to ensure data quality for information analysis, including removing contaminated reads from joints (the number of nucleobase > 5 bp), removing low quality reads (Phred Quality Score ≤ 19 account for 15% of total nucleobase) and reads with N ratio greater than 5%. The data and quality of the clean reads were calculated, including Q30 and GC-content (Supplementary Table S1).
Transcriptome assembly and gene annotation. Clean reads were assembled using Trinity software (v2.4.0) 42 to obtain full-length transcripts (Supplementary Table S2 Table S4). All DEGs were mapped to GO terms to characterize main biological functions (False discovery rate, FDR < 0.05). Kyoto Encyclopedia of Genes and Genomes (KEGG, http:// www. genome. jp/ kegg) pathway was used to analyze DEGs and identify the enriched metabolic pathways. GO and KEGG enrichments were analyzed using Hyper geometric test. mRNA library validation. Twelve highly expressed DEGs were randomly selected for validation (including three up-regulated and nine down-regulated DEGs). The 18 S gene was used as the reference gene and the annealing temperature was set at 60 °C 46,47 . The DEGs and 18 S were used for the validation of RNA-seq data by real-time PCR ( Table 2). The PrimeScript™ RT reagent Kit (TaKaRa, Japan) was used for reverse transcription of the RNA of sea urchins at three different flow velocities according to the instruction and system. Fluorescence quantitative analysis was performed in Light Cycler 96 and the fluorescence quantitative kit was SYBR ® Premix Ex TAq™ (TaKaRa, Japan). The qRT-PCR reaction system had a total volume of 20 µL, including 2 µL of cDNA template, 10 µL of 2 × SYBR Green Master mix (TaKaRa, Japan), 0.8 µL of each primer and 6.4 µL of PCR-grade water. The running program was set as follows: 95 °C for 30 s and followed by 40 cycles. The annealing and elongation were 95 °C for 5 s and 60 °C for 32 s, respectively. PCR melting curve analysis was conducted to confirm single PCR products at the end of the reaction. The fluorescence quantitative results were calculated using the 2 −ΔΔCT algorithm 48 . The results of each sample were tested three times.
Statistical analysis. Normal distribution and homogeneity of variance were analyzed using Kolmogorov-Smirnov test and Levene test, respectively. Gene expressions were analyzed using independent-samples t-test. All statistical analyses were performed by using SPSS 22.0 (Inc., Chicago, IL, USA). All data were expressed as mean values ± standard deviation (mean ± SD) and evaluated the significance at the level of P < 0.05. Ethical approval. All applicable international, national, and/or institutional guidelines for the care and use of animals were followed by the authors.

Data availability
The raw data were submitted to the NCBI database (Accession Number: PRJNA825640).